function [] = testLGC()

    %% initialization:
    iSeqLength = 10000;
    x_i = zeros(1, iSeqLength); % create empty sequence for the pseudorandom numbers ...
    u_i = x_i; % initialize the sequence for the pseudorandom rationals ...

    % constant parameters ...
    m = 2^35; % modulus-value, 0 < m.
    a = 2^18; % multiplier, 0 < a < m.
    x_0 = 314159265; % initial start value ("seed"), 0 <= x_0 < m. 
        
    x_i(1) = x_0; % initialization of the sequence ...

    %% generate pseudorandom numbers between the interval (0, 1) with the
    %% help of the "classic" Linear Congruential Generator (LCG):
    for i = 1:iSeqLength-1
        x_i(i+1) = mod( ( (a+1)*x_i(i) + 1 ), m );
    end
    % map the generated values into pseudorandom rationals u_i in the
    % intervall (0, 1) ...
    u_i = x_i./m;
    
    %% grafical outputs of the value distributions:
    
    cGrey = [0.4, 0.4, 0.4];
    cDarkGrey = [0.3, 0.3, 0.3];
    cLightGrey = [0.6, 0.6, 0.6];
    iStdSize = 12;
    iTitSize = 14;

    figure(1); % 2D-plot ...
    clf;

    hold on;
    for i = 1:iSeqLength-2
        plot(u_i(i), u_i(i+1), '.', 'Color', cDarkGrey);
    end
    xlbl = xlabel('$u_i$');
    ylbl = ylabel('$u_{i+1}$');
    set(xlbl, 'Interpreter', 'latex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(ylbl, 'Interpreter', 'latex', 'FontName', 'Times', 'FontSize', iStdSize);
    title('2D:', 'FontName', 'Times', 'FontSize', iTitSize);
    
    axis square;    
    grid on;
    hold off;

    figure(2); % 3D-plot ...
    clf;

    hold on;
    for i = 1:iSeqLength-2
        plot3(u_i(i), u_i(i+1), u_i(i+2), '.', 'Color', cGrey);
    end
    xlbl = xlabel('$u_i$');
    ylbl = ylabel('$u_{i+1}$');
    zlbl = zlabel('$u_{i+2}$');
    set(xlbl, 'Interpreter', 'latex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(ylbl, 'Interpreter', 'latex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(zlbl, 'Interpreter', 'latex', 'FontName', 'Times', 'FontSize', iStdSize);
    title('3D:', 'FontName', 'Times', 'FontSize', iTitSize);

    axis square;
    grid on;
    view(-25, 45);
    hold off;
    
    %% Pearson's Chi-square test:
    % to verify the the distribution property of a given sample set ...
        
    % parameter initialization:
    iNObs = 17 % number of observation classes ...
    %iNObs = 5
    O_i = zeros(1, iNObs); % init. observation class list ...
    E_i = O_i; % init. of the expected frequencies array ...
    
    alpha = 0.05   % probability that Chi^2 exceeds
                   % critical value c = P(Chi^2 >= c).
                    
    % critical value c for the alpha = 0.05 and the degree of 
    % freedom df = 4 (iNObs - 1).
    % compare given table in: Schaum's Outline, Probability, Second Edition, p. 284.
    c = 9.49;

    % count the number of values of each observation class (frequency) ...
    [O_i, xout] = hist(u_i, iNObs);
    
    % optional ...
    % put out a small difference than the hist-function of matlab.
    % hist:                 (2016, 1969, 2015, 2007, 1993)
    % count func. below:    (2015, 1970, 2015, 2007, 1993)
    %
    % calculate the slot size of each class ...
    %iStepSize = 1.0/iNObs
    %
    %iL = 0;
    %iR = iStepSize;
    %for i = 1:iNObs
    %   for j = 1:iSeqLength
    %      if ( (u_i(j) >= iL) && (u_i(j) < iR) )
    %          O_i(i) = O_i(i) + 1; 
    %      end
    %   end
    %   % next class ...
    %   iL = iR;
    %   iR = iR + iStepSize;
    %end    
   
    % plot the number of values of each class (histogramm) ...
    figure(3);
    clf;
    bar(xout, O_i);
    
    h = findobj(gca,'Type','patch');
    set(h,'FaceColor', cLightGrey,'EdgeColor','k')
    
    xlbl = xlabel('\it u_i');
    ylbl = ylabel('\it Frecuencia');
    set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
    
    % nested function to verify the null hypothesis of the 
    % goodness-of-fit tests (Chi-square and KS):
    function checkNullHypothesis(stat, c)
        if (stat >= c)
            fprintf('The test does not fit!\nRejecting the null hypothesis!\n');
            fprintf('(%f >= %f)\n\n', stat, c);
        else
            fprintf('The value does not exceed the critical value c (%f < %f).\n', stat, c);
            fprintf('Null hypothesis accepted!\n\n');
        end
    end


    %% Chi-square test against the normal distribution (PDF):
    % note: the result never fits! For the given alpha = 0.05, it exceeds
    % massively all the critical values c of the table in the book
    % [Schaum's Outline, Probability, Second Edition, p. 284].
    disp('Chi-square test against the normal distribution (PDF):');
    
    % calculate the arithmetic average of X (mu = E(X) = sum(x_i*f(x_i), i=1, n), f(x_i) = 1/n ) 
    % by assuming a discrete uniform distribution ...
    mu = mean(O_i)
    %mu = sum(O_i)/iNObs; % optional ...
    
    % calculate the sample variance S^2 ...
    S = std(O_i, 1)
    %S2 = S^2;
    %S2 = sum( (O_i - mu).^2 ) * 1.0/iNObs; % optional ...
    %S2 = sum(O_i.^2)/iNObs - mu^2; % optional ...
   
    % calculate expected (theoretical) frequency, asserted by the
    % null hypothesis ...
	E_i = iNObs * normpdf(O_i, mu, S)
    %E_i = iNObs * ( 1.0/sqrt(2*pi*S2) * exp(-(O_i - mu).^2 / (2.0*S2)) ); % optional ...
    
    % plot the expected frequency ...
    figure(4);
    clf;
    
    plot(   O_i, E_i, 'ko', 'LineWidth', 1, 'MarkerEdgeColor', 'k', ...
            'MarkerFaceColor', cLightGrey, 'MarkerSize', 5  );
    xlbl = xlabel('\it O_i');
    ylbl = ylabel('\it E_{i}');
    set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
    %title( sprintf( 'Distribuci%cn de probabilidad del valores esperado te%crico', char(162), char(162)), ...
    %        'FontName', 'Times', 'FontSize', iTitSize  );
    axis normal;
    grid on;
    
    chi2_0 = sum( ((O_i - E_i).^2) ./ E_i)

    % calculate the degree of liberty df ...
    df = iNObs - 1 - 2    % + (-2) because other statistics are
                          % additionally used (mean & variance) ...
    c = getCriticalValue_Chi2(alpha, df)
    % check the null hypothesis ...
    checkNullHypothesis(chi2_0, c);
    

    %% Chi-square test against the uniform distribution:
    disp('Chi-square test against the uniform distribution:');
    
    % calculate expected (theoretical) frequency, asserted by the 
    % null hypothesis ...
    E_i(1, :) = iSeqLength/iNObs
        
    % calculate the value of the test statistic (Chi^2) against the 
    % standard uniform distribution ...
    chi2_0 = sum( ((O_i - E_i).^2) ./ E_i)
            
    % calculate the degree of liberty df ...
    df = iNObs - 1
    % get the critical value of the Chi-square distribution table ...
    c = getCriticalValue_Chi2(alpha, df)
    % finally verify the null hypothesis ...
    checkNullHypothesis(chi2_0, c);
    
%     % alternative method with matlab function (result does not match) ...
%     edges = linspace(0, 1, iNObs+1)
%     [h,p,st] = chi2gof(O_i, 'edges', edges, 'expected', E_i)

        
    %% Kolmogorov-Smirnov test:
    disp('Kolmogorov-Smirnov test:');
    
    n = 10; % number of observation classes ...
    
    % initialization:
    F_n = zeros(1, n); % value array for the empirical distribution function F_n(x) ...
    F_x = F_n; % init. value array for the CDF F(x) of the real sample values ...
         
    % calculate the values of the empirical distribution function F_n(x)
    % for n (independent identically distributed) observations X_i
    % s.t.: F_n(x) = 1/n * sum(I(X_i <= x), i=1, n), where
    %       I(X_i <= x) is the indicator function.
    %
    %                     |0       x < X_1
    %           F_n(x) = -|i/n     X_i <= x < X_{i+1}
    %                     |1       x >= X_n
    %
    % it gives the relative frequency with which the observed values are <= x. 
    % [counting up the number of elements in the sample (cumulatively) and divide by n.]
    F_n(1, :) = [1:n]./n    % in this case the cumul. distribution function
                            % is linear (because we assume a uniform distrib.).
    
    % create the cumulative distribution function F(x) = x
    % (build the CDF from the real sample values).
    [F_x, x_val] = hist(u_i, n)
    % optional ...
    % create an interval vector with 0 <= x <= 1
%     x_val = linspace(0, 1, n+1) % vector with constant intervals ...
%     for i = 1:n
%         for j = 1:iSeqLength
%             %if ( (u_i(j) <= x_val(i)) && (x_val(i) < u_i(j+1)) )
%             if ( (x_val(i) < u_i(j)) && (u_i(j) <= x_val(i+1)) )
%                 F_x(i) = F_x(i) + 1;
%             end
%         end
%     end
      
    % plot a histogramm of the cumulative values ...
    figure(5);
    clf;
    bar(x_val(1:n), F_x);
    
    h = findobj(gca,'Type','patch');
    set(h,'FaceColor', cLightGrey, 'EdgeColor', 'k');
    
    xlbl = xlabel('\it x');
    ylbl = ylabel('\it F(x)');
    set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
    set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);

    F_x = F_x./iSeqLength; % normalizing ...
    % sort the values in ascending order ...
    F_x = sort(F_x)
         
    % calculate KS-statistics (difference of the distribution functions).
    % searching for the largest difference between F_n and F:
    D_plus = max(F_n - F_x) % pos. difference ...
    D_minus = max(F_x - [0:n-1]./n) % neg. difference ...
    
    D = max(D_plus, D_minus) % max. difference value ...
    
    % get the critical value ...
    c = getCriticalValue_KS(alpha, n)
    % verify the null hypothesis ...
    checkNullHypothesis(D, c);
    
    % plot the cumulative distribution functions ...
    figure(6);
    clf;
       
%     xc_a = zeros(1, n);
%     for i = 1:n
%        xc_a(i) = sum(F_x(1:i)); 
%     end

    x = linspace(0, 1, n+1);
%     f = stairs([x_val(1:n)], [xc_a], '-');
    f = cdfplot(F_x);
    hold on
    g = plot(x, [0, F_n], 'k-');
    
    set(f, 'LineWidth', 1, 'Color', cDarkGrey);
    set(g, 'LineWidth', 2);
    
    l = legend( [f g], 'CDF $F(x)$', 'Empirico-lineal $F_n(x)$' );
    set( l, 'Interpreter', 'latex', 'Location', 'NorthWest', ...
         'FontName', 'Times', 'FontSize', 9 );
    % correct the box size of the legend (make it larger) ...
    pos = get(l, 'position');
    pos(3) = 1.2*pos(3);
    set(l, 'position', pos);
    
%    axis([0 1 0 1]);
   axis([0.096 0.102 0 1]);
   hold off;
   
   % matlab function for SM-test
   % (to verify the result of manually calculated version) ...
   disp('using the given function of Matlab (kstest):');
   [h, p, ksstat, cv] = kstest(F_x, [F_x', F_x'], alpha)
   % check null hypothesis ...
   checkNullHypothesis(ksstat, cv); % get the same result!!! :-)
     
%    [f, x, flo, fup] = ecdf(F_x)%, 'frequency', x_val)    
%    figure(7);
%    clf;
%    
%    stairs(f, x, 'r')
%    %hold on;
%    %stairs(f, flo, 'g');
%    %stairs(f, fup, 'b');
%    %hold off;
    
end

